Seismic Data Processing

ABSTRACT

The invention includes a method for reducing noise in migration of seismic data, particularly advantageous for imaging by simultaneous encoded source reverse-time migration (SS-RTM). One example embodiment includes the steps of obtaining a plurality of initial subsurface images; decomposing each of the initial subsurface images into components; identifying a set of components comprising one of (i) components having at least one substantially similar characteristic across the plurality of initial subsurface images, and (ii) components having substantially dissimilar characteristics across the plurality of initial subsurface images; and generating an enhanced subsurface image using the identified set of components. For SS-RTM, each of the initial subsurface images is generated by migrating several sources simultaneously using a unique random set of encoding functions. Another embodiment of the invention uses SS-RTM for velocity model building.

CROSS-REFERENCE TO RELATED APPLICATION

This application is a continuation-in-part of U.S. application Ser. No.13/430,380, filed Mar. 26, 2012, entitled SEISMIC DATA PROCESSING andclaims the benefit of U.S. Provisional Patent Application 61/479,199,filed Apr. 26, 2011, entitled SEISMIC DATA PROCESSING, the entirety ofwhich is incorporated by reference herein.

FIELD OF THE INVENTION

The invention relates generally to the field of geophysical prospecting,and more particularly to seismic data processing. Specifically, theinvention is a method for identifying and/or reducing signal noiseintroduced by various conventional inversion techniques.

BACKGROUND OF THE INVENTION

Geophysical inversion attempts to find a model of subsurface propertiesthat optimally explains observed data and satisfies geological andgeophysical constraints. There are a large number of well known methodsof geophysical inversion. These well known methods fall into one of twocategories, iterative inversion and non-iterative inversion. Thefollowing are definitions of what is commonly meant by each of the twocategories:

-   -   a. Non-iterative inversion—inversion that is accomplished by        assuming some simple background model and updating the model        based on the input data. This method does not use the updated        model as input to another step of inversion. For the case of        seismic data these methods are commonly referred to as imaging,        migration, diffraction tomography or Born inversion.    -   b. Iterative inversion—inversion involving repetitious        improvement of the subsurface properties model such that a model        is found that satisfactorily explains the observed data. If the        inversion converges, then the final model will better explain        the observed data and will more closely approximate the actual        subsurface properties. Iterative inversion usually produces a        more accurate model than non-iterative inversion, but is much        more expensive to compute.

In general, inversion is beneficial in correcting observed seismic dataso that reflections are plotted at a true representation of theirsubsurface locations (Stolt, 1978; Claerbout, 1985). The need to correct(i.e., invert) observed seismic data arises, for example, becausereflections from dipping interfaces are observed, and thereforerecorded, at surface positions that are not directly above thesubsurface locations where the reflections actually occur. Also,isolated point-like discontinuities in the subsurface (i.e., pointscatterers) generate seismic events (e.g., diffractions) recorded over alarge range of receivers. Such diffractions can make the properinterpretation of seismic data confusing. Furthermore, seismic velocityvariations can also cause a distorted view of the subsurface geology(Yilmaz, 1987). It is only after inversion that the structures andgeometric configurations observed in seismic recordings can be thoughtof as an accurate depiction of the geologic layers that gave rise to theseismic reflections.

Because of the increased complexity of iterative inversion, it isgenerally desirable to use a non-iterative form of inversion (i.e.,imaging, migration, etc.) when possible. However, as industry exploresmore complex geographic areas, traditional imaging and interpretationmethods fail to provide subsurface images having the quality (e.g.,accuracy) desired in making decisions on exploration and production. Forexample, wave-equation migration algorithms are based upon the one-waywave-equation approach. The one-way wave equation assumes that wavespropagate in only one primary direction, either down into the subsurfaceor up from the subsurface. Because of the one-directional nature ofpropagation, imaging steeply dipping reflectors is difficult.

Advanced imaging tools that use full physics of wave-propagation, suchas reverse time migration (i.e., RTM), generally provide better imagesof the subsurface. Such approaches use solutions of the two-way waveequation. Migration techniques that use the two-way wave equationgenerally provide a more accurate result because waves propagating inall directions are handled equally well, and wave amplitudes areproperly modeled since no approximations are used in the algorithm.However, there is a cost associated with conventional two-way waveequation techniques. Specifically, the full physics of propagating wavesin a complex geological setting where the medium velocity is complicatedis very computationally intensive. The computational demand is furtheraccentuated when there are many (e.g., thousands) shot records to bemigrated in a three-dimensional (3D) setting and/or when high frequencydata is obtained in an effort to increase the resolution of thesubsurface images. It is desirable, therefore, to have a system and/ormethod for increasing the computational efficiency (i.e., reducing thecomputational cost) of two-way wave equations based techniques, such asRTM.

SUMMARY OF THE INVENTION

The present disclosure relates generally to the field of geophysicalprospecting, and more particularly to seismic data processing. Oneexemplary embodiment includes the steps of obtaining a plurality ofinitial subsurface images; decomposing each of the initial subsurfaceimages into components; identifying a set of components comprising oneof (i) components having at least one substantially similarcharacteristic across the plurality of initial subsurface images, and(ii) components having substantially dissimilar characteristics acrossthe plurality of initial subsurface images; and generating an enhancedsubsurface image using the identified set of components. Each of theinitial subsurface images is generated using a unique random set ofencoding functions.

Another exemplary embodiment includes the steps of obtaining a pluralityof simultaneous-source reverse-time migration (SS-RTM) subsurfaceimages; decomposing each of the plurality of SS-RTM subsurface imagesinto curvelet coefficients; averaging the curvelet coefficients togenerate a preliminary signal curvelet coefficient estimate; computing avariance of a subset of the curvelet coefficients to determine a noiselevel in the preliminary signal curvelet coefficient estimate;attenuating noise in the curvelet coefficients using the determinednoise level and the preliminary signal curvelet coefficient estimate togenerate attenuated curvelet coefficients; and performing an inversecurvelet transform on the attenuated curvelet coefficients to generatean enhanced SS-RTM image. Each of the SS-RTM subsurface images isgenerated using a unique random unit-magnitude set of encodingfunctions.

Yet another exemplary embodiment includes computing a gradient of a costfunction associated with seismic data by generating a plurality ofgradients of objective functions computed using unique sets of randomreciprocal encoding functions; decomposing each of the gradients intocomponents; identifying a set of components having at least onesubstantially similar characteristic across the plurality of gradients;and generating an enhanced gradient using the identified set ofcomponents.

Another example embodiment includes using SS-RTM for velocity modelbuilding. An advantage of using SS-RTM is that it speeds up the velocitymodel building process. In a typical migration algorithm, the forwardand the backward wave field are cross-correlated to form an image. Ifthe velocity model is accurate, then the cross-correlation at zero-lagcross-correlation will focus the energy constructively to form theimage. If the velocity is incorrect then the focusing is poor. Intraditional practice, image gathers are generated and flatness of thesegathers is examined as a metric to qualify the accuracy of the velocitymodel. However, the process of generating image gathers from wave-basedmigration (such as RTM) is an expensive process particularly whenindividual shots are used in the analysis. Typically in a 3D survey onehas several thousands of shots. This embodiment of the present inventivemethod offers a cheaper alternative that can aid in judging the qualityof the velocity model.

BRIEF DESCRIPTION OF THE DRAWINGS

The present invention will hereinafter be described in conjunction withthe following figures, wherein like numerals denote like elements, and

FIG. 1 is a flow diagram of a method for iteratively performing SS-RTMin accordance with an embodiment of the present invention;

FIG. 2 is a flow diagram of a method for generating an enhanced image inaccordance with an embodiment of the present invention;

FIGS. 3A-4D are images corresponding to an exemplary simulatedapplication of one or more embodiments of the present invention;

FIGS. 4A-4B illustrate the distribution in different domains of thenoise corresponding to the exemplary simulated application of FIG. 3.

FIG. 5 shows sequential-shot RTM images of seismic data binned intodifferent offset bins;

FIG. 6 shows SS-RTM of the same seismic data as in FIG. 5, binned intodifferent offset bins; and

FIG. 7 shows SS-RTM images for the four offset bins with an incorrect(slow) velocity model.

Due to patent law restrictions on use of color, FIGS. 4A-4B are blackand white reproductions of color drawings.

DETAILED DESCRIPTION Definitions

Various terms as used herein are defined below. To the extent a termused in a claim is not defined below, it should be given the definitionpersons in the pertinent art have given that term.

As used herein, the “a” or “an” entity refers to one or more of thatentity. As such, the terms “a” (or “an”), “one or more”, and “at leastone” can be used interchangeably herein unless a limit is specificallystated.

As used herein, the terms “comprising,” “comprises,” “comprised,” and“comprise” are open-ended transition terms used to transition from asubject recited before the term to one or more elements recited afterthe term, where the element or elements listed after the transition termare not necessarily the only elements that make up of the subject.

As used herein, the terms “containing,” “contains,” and “contain” havethe same open-ended meaning as “comprising,” “comprises,” and“comprise.”

As used herein, the term “exemplary” means “serving as an example,instance, or illustration.” Any embodiment described herein as“exemplary” is not necessarily to be construed as preferred oradvantageous over other embodiments.

As used herein, the terms “having,” “has,” and “have” have the sameopen-ended meaning as “comprising,” “comprises,” and “comprise.”

As used herein, the terms “including,” “includes,” and “include” havethe same open-ended meaning as “comprising,” “comprises,” and“comprise.”

As used herein, the term “substantially” refers to the complete ornearly complete extent or degree of an action, characteristic, property,state, structure, item, or result. For example, an object that is“substantially” enclosed would mean that the object is either completelyenclosed or nearly completely enclosed. The exact allowable degree ofdeviation from absolute completeness may in some cases depend on thespecific context. However, generally speaking the nearness of completionwill be so as to have the same overall result as if absolute and totalcompletion were obtained. The use of “substantially” is equallyapplicable when used in a negative connotation to refer to the completeor near complete lack of an action, characteristic, property, state,structure, item, or result.

Description

For illustrative purposes, the present invention is described below inconnection with reverse time migration and, more specifically,simultaneous source reverse time migration (i.e., SS-RTM). However, itwill be apparent to one skilled in the art having the benefit of thisdisclosure that the present technique may be used in connection with anytwo-way wave inversion technique.

The significant computational cost of RTM imaging is driven by thewavefield extrapolation procedure. This procedure comprises accurateforward propagation of the source wavefield as well as backwardpropagation of the recorded receiver wavefield. These propagatedwavefields are then combined by an imaging condition (e.g.,cross-correlation) to produce a corresponding subsurface image(Whitmore, 1983; Baysal et al., 1983). In traditional depth-imaging ofseismic data, individual shot gathers (i.e., source and receiverwavefield pairs corresponding to source and receiver hardware pairs) areimaged separately; that is, for each shot, the source and receiverwavefields are extrapolated and then combined to produce the subsurfaceimage. The complete RTM image is constructed by combining the subsurfaceimages generated from each shot. For such an approach, the cost ofimaging is directly proportional the number of shots. Thus the cost ofRTM increases significantly as the number of shots increases.

Simultaneous source migration (Romero et al., 2000) is one techniquethat may be used to reduce the computational cost of RTM. Insimultaneous source imaging, several shots are encoded using a randomencoding function to form a super-shot, and this super-shot is thenmigrated. The computational cost of migrating a super-shot is similar tothe cost of migrating a single shot. Thus the use of simultaneous sourcemigration may significantly reduce the cost of RTM. Unfortunately, theuse of simultaneous source migration generally introduces undesirablenoise into the resulting subsurface image. Most notably, the imagesuffers from interference between the different sources. Consequently,the resulting SS-RTM image suffers from low signal to noise ratio. Thepresent inventive method introduces a denoising strategy to address thisproblem.

The basic mathematical framework of simultaneous source migration can bedescribed as explained below. Let S(i) denote each source wavefield, ofwhich there is generally many, and R(i,j) denote the receiver wavefieldactivated by a given S(i). The interference arises because thecross-correlation between the forward propagated S(i) andback-propagated R(k,j), for k≠i, is not zero. More specifically:

Effective SuperShot=Σ_(i) a(i)S(i)   [EQ 1]

Effective SuperShot receiver measurement at j-th receiver=Σ_(i)a(i)R(i,j)   [EQ2]

In conventional SS-RTM, the Effective SuperShot source would be forwardpropagated, and the effective SuperShot receiver measurement would beback-propagated. The results of the forward and back propagation wouldthen be crosscorrelated (denoted by {circle around (×)}). Let F denoteforward propagation and B denote back propagation operators.

$\begin{matrix}\begin{matrix}{{Image} = {{F\left\lbrack {\sum\limits_{i}{{a(i)}{S(i)}}} \right\rbrack} \otimes {B\left\lbrack {\sum{{a(k)}{R\left( {k,j} \right)}}} \right\rbrack}}} \\{= {\left( {\sum\limits_{i}{{{a(i)}}^{2}{{F\left\lbrack {S(i)} \right\rbrack} \otimes {B\left\lbrack {R\left( {i,j} \right)} \right\rbrack}}}} \right) +}} \\{\left( {\sum\limits_{i,k}{{a(i)}a*(k){{F\left\lbrack {S(i)} \right\rbrack} \otimes {B\left\lbrack {R\left( {k,j} \right)} \right\rbrack}}}} \right)}\end{matrix} & \lbrack{EQ3}\rbrack\end{matrix}$

where * denotes complex conjugation.

-   Note that:

Signal (1st term in Image)=F[S(i)]

B[R(i,j)] and   [EQ4]

Interference Noise (2nd term in Image)=F[S(i)]

B[R(k,j)], with k≠i   [EQ5]

It may be beneficial, therefore, to distinguish between, and thensubsequently separate, the signal and noise components from the SS-RTMimage. It may be noted that instead of measured receiver wavefields andsource wavefields; one can use preconditioned/pre-processed receiver andsource wavefields to construct the SS-RTM image.

One conventional technique that attempts to reduce the noise componentis to use random source encoding (Romero et al., 2000; Perrone and Sava,2009, 2010). Typically several runs with different encodings are carriedout and the images from those runs are then stacked; see Ober et al,U.S. Pat. No. 6,021,094, and Romero et al, 2000. Such a method mayproduce an image with less artifacts. However, such an approach reliessolely on averaging to perform noise attenuation and does not fullyleverage the diversity of information in the available encoded SS-RTMimages.

Another conventional technique that may reduce cross-talk artifacts isto formulate the Simultaneous Source imaging problem as a least-squaresmigration problem (Verschuur and Berkhout, 2009; Tang and Biondi, 2009;Dai and Schuster, 2009). Perrone and Sava (2009; 2010) investigateddifferent encoding schemes as well as mixed encoding schemes such as acombination of linear and random encoding to reduce cross-talkartifacts. Godwin and Sava (2010) use singular value decomposition todetermine an encoding scheme that attempts to balance computational costsavings with the number of cross-talk artifacts. The Godwin-Sava methodinherently assumes that the images generated using each encoding isstacked. As such the Godwin-Sava method does not permit the estimationof the signal components that are not imaged in the final result withoutconducting additional SS-RTM runs and stacking the resulting images.

It may be appreciated from the above discussion, then, that a needexists for a better system and/or method for reducing noise component(i.e., cross-talk artifacts). The present invention provides such asystem and/or method. In contrast to conventional techniques whichgenerally use the same encoding for each source and receiver pair (seeEQ1 and EQ2), the present invention may, in at least one embodiment,optionally use a first set of encodings for the sources and a second setof encodings for the receivers. That is:

Effective SuperShot=Σ_(i) a(i)S(i), and   [EQ6]

Effective SuperShot receiver measurement at j-th receiver=Σ_(i)b(i)R(i,j)   [EQ7]

Consequently, the effective SS-RTM image may be given by:

$\begin{matrix}\begin{matrix}{{Image} = {{F\left\lbrack {\sum\limits_{i}{{a(i)}{S(i)}}} \right\rbrack} \otimes {B\left\lbrack {\sum{{b(k)}{R\left( {k,j} \right)}}} \right\rbrack}}} \\{= {\left( {\sum\limits_{i}{{a(i)}b*(i){{F\left\lbrack {S(i)} \right\rbrack} \otimes {B\left\lbrack {R\left( {i,j} \right)} \right\rbrack}}}} \right) +}} \\{\left( {\sum\limits_{i,k}{{a(i)}b*(k){{F\left\lbrack {S(i)} \right\rbrack} \otimes {B\left\lbrack {R\left( {k,j} \right)} \right\rbrack}}}} \right)}\end{matrix} & \lbrack{EQ8}\rbrack\end{matrix}$

In at least one other embodiment, however, the present invention may usethe same and/or equivalent encoding for each source and receiver pair.

In general, the present invention comprises the following steps toimprove signal and interference noise separation in SS-RTM images.First, SS-RTM is performed a plurality of times to generate a pluralityof subsurface images. Different “reciprocal” sets of source and receiverrandom encoding functions are used during each application of SS-RTM tothe set of shot gathers. It may be understood that a pair of randomencoding functions, such as a_(n)(i)and b_(n)(i), are “reciprocal” whena_(n)(i)=1/b_(n)*(i) or a_(n)(i)=b_(n)(i)=0, where n denotes the SS-RTMrun index. Second, the resulting signal and noise characteristics areused along with the plurality of resulting SS-RTM images to generate anenhanced image with reduced interference noise. With regard to thesecond step, and as previously discussed, one conventional strategy toobtain an enhanced image is to average (i.e., stack) a plurality ofgenerated subsurface images; see Ober et al, U.S. Pat. No. 6,021,094,and Romero et al, 2000. During stacking, the signal component addsconstructively (i.e., stacks in) and the noise component gets attenuated(i.e., stacks out). The stacked image can be mathematically expressed as

$\begin{matrix}{{Stack} = {\frac{\sum\limits_{n = 1}^{N}{Image}_{n}}{N} = {{Signal} + \frac{\sum\limits_{n = 1}^{N}{Noise}_{n}}{N}}}} & \lbrack{EQ9}\rbrack\end{matrix}$

The signal-to-noise-ratio (i.e., SNR) of the signal and noise events maybe predicted as a function of the number of runs. The noise energy wouldbe controlled by the average of a string of random phase (such as ±1)numbers. It is possible to show that the all interference noise termswould cancel out after stacking only when N_(s) SS-RTM runs areconducted, with N_(s) denoting the number of shots, and the encodingschosen to be orthogonal to each other. Note that if N_(s) SS-RTM runsare conducted, then no efficiency is gained. When substantially fewerthan N_(s) SS-RTM runs (i.e., N<<N_(s)) are conducted, then simplestacking of the SS-RTM images does not adequately attenuate theinterference noise. It is desirable, therefore, to have a system and/ormethod for generating an enhanced image that does not exhibit thedeficiencies of the stacking method.

With reference now to FIG. 1, a flow diagram of a method 100 foriteratively performing SS-RTM is provided in accordance with at leastone embodiment of the present invention. The method 100 may beadvantageously implemented in connection with any appropriate methodand/or system, such as a computer implemented system, to meet the designcriteria of a particular application. The method 100 generally includesa plurality of blocks or steps that may be performed serially. As willbe appreciated by one of ordinary skill in the art, the order of thesteps shown in FIG. 1 is exemplary and the order of one or more stepsmay be modified within the spirit and scope of the present invention.Additionally, the steps of the method 100 may be performed in at leastone non-serial (or non-sequential) order, and one or more steps may beomitted to meet the design criteria of a particular application. Block102 represents an entry point into the method 100.

At block 104, a Set_(n) of Random Encoding Functions is selected. TheSet_(n) comprises Random Encoding Functions A_(n) and B_(n), where B_(n)is generally reciprocal to A_(n). It may be noted, however, that in atleast one embodiment A_(n) may be identical to or equivalent to B_(n.)

At block 106, the Forward (i.e., source) Wave Components of a set ofshot gathers are encoded using the Random Encoding Function A_(n) togenerate an Encoded Source Super-Shot Wave Component.

At block 108, the Backward (i.e., received) Wave Components of the setof shot gathers are encoded using Random Encoding Function B_(n) togenerate an Encoded Receiver Super-Shot Wave Component.

At block 110, the Encoded Source Super-Shot Wave Component is forwardpropagated to generate a Forward Propagated Wave Component.

At block 112, the Encoded Receiver Super-Shot Wave Component is backpropagated to generate a Back Propagated Wave Component.

At block 114, an imaging condition is applied to the Forward PropagatedWave Component and the Back Propagated Wave Component to generate aSubsurface Image (SI_(n)).

The steps represented by blocks 104-114 may be performed for anyappropriate number of iterations to satisfy the design criteria of aparticular application. Decision block 116 represents the determinationof whether or not additional iterations are desirable. If additionaliterations are desirable, then the method 100 generally returns toblock/step 104. Each Set_(n) selected at block 104 is unique in thateach set includes at least one member that makes the set, as a whole,different from any previously selected set. For example, a Set_(n) mayinclude Reciprocal Random Encoding Functions X and Y and Set₂ mayinclude Reciprocal Random Encoding Functions X and Z. Set₁ and Set₂ areunique because the set XY is not identical to the set XZ. It may beappreciated that subscript n represents an indexing digit that isgenerally modified upon each iteration.

If the method 100 does not return to block 104 then the method 100generally falls through to block 118. Block 118 represents an exit fromthe method 100.

As illustrated in FIG. 1, then, SS-RTM is generally performed repeatedly(i.e., more than once) on a given seismic data set (e.g., shot gatherset) using reciprocal sets of source and receiver random encodingfunctions a_(n)(i)'s and b_(n)(i)'s respectively. For such reciprocalencodings, the n-th SS-RTM image would be described by:

$\begin{matrix}\begin{matrix}{{Image}_{n} = {{F\left\lbrack {\sum\limits_{i}{{a_{n}(i)}{S(i)}}} \right\rbrack} \otimes {B\left\lbrack {\sum{{b_{n}(k)}{R\left( {k,j} \right)}}} \right\rbrack}}} \\{= {\left( {\sum\limits_{i}{{a_{n}(i)}b_{n}*(i){{F\left\lbrack {S(i)} \right\rbrack} \otimes {B\left\lbrack {R\left( {i,j} \right)} \right\rbrack}}}} \right) +}} \\{\left( {\sum\limits_{i,k}{{a_{n}(i)}b_{n}*(k){{F\left\lbrack {S(i)} \right\rbrack} \otimes {B\left\lbrack {R\left( {k,j} \right)} \right\rbrack}}}} \right)} \\{= {\left( {\sum\limits_{i}{{F\left\lbrack {S(i)} \right\rbrack} \otimes {B\left\lbrack {R\left( {i,j} \right)} \right\rbrack}}} \right) +}} \\{= {\left( {\sum\limits_{i}{{F\left\lbrack {S(i)} \right\rbrack} \otimes {B\left\lbrack {R\left( {i,j} \right)} \right\rbrack}}} \right) +}} \\{{\left( {\sum\limits_{i,k}{{a_{n}(i)}b_{n}*(k){{F\left\lbrack {S(i)} \right\rbrack} \otimes {B\left\lbrack {R\left( {k,j} \right)} \right\rbrack}}}} \right),}} \\{= {{Signal} + {Noise}_{n}}}\end{matrix} & \lbrack{EQ10}\rbrack\end{matrix}$

The use of reciprocal encodings insures that each SS-RTM image iscomprised of both the desired signal (with unit magnitudes in eachImage) and a noise that varies across the different SS-RTM images.

In at least one embodiment, unit-magnitude complex numbers may be usedas a_(n)(i)'s. In such an embodiment a_(n)(i)=b_(n)(i),|a_(n)(i)|=|b_(n)(i)|=1, and the n-th SS-RTM image would be describedby:

$\begin{matrix}\begin{matrix}{{Image}_{n} = {{F\left\lbrack {\sum\limits_{i}{{a_{n}(i)}{S(i)}}} \right\rbrack} \otimes {B\left\lbrack {\sum{{a_{n}(k)}{R\left( {k,j} \right)}}} \right\rbrack}}} \\{= {\left( {\sum\limits_{i}{{{a_{n}(i)}}^{2}{{F\left\lbrack {S(i)} \right\rbrack} \otimes {B\left\lbrack {R\left( {i,j} \right)} \right\rbrack}}}} \right) +}} \\{{\left( {\sum\limits_{i,k}{{a_{n}(i)}a_{n}*(k){{F\left\lbrack {S(i)} \right\rbrack} \otimes {B\left\lbrack {R\left( {k,j} \right)} \right\rbrack}}}} \right),}} \\{= {\left( {\sum\limits_{i}{{f\left\lbrack {S(i)} \right\rbrack} \otimes {B\left\lbrack {R\left( {i,j} \right)} \right\rbrack}}} \right) +}} \\{{\left( {\sum\limits_{i,k}{{\exp \left( {j\; {\theta_{n}\left( {i,k} \right)}} \right)}{{F\left\lbrack {S(i)} \right\rbrack} \otimes {B\left\lbrack {R\left( {k,j} \right)} \right\rbrack}}}} \right),}} \\{= {{Signal} + {Noise}_{n}}}\end{matrix} & \lbrack{EQ11}\rbrack\end{matrix}$

where θ_(n)(i,k) is the phase difference between a_(n)(i) and a_(n)*(k).Once again, the use of unit-magnitude encodings generally provides asignal component that does not vary and a noise component that varies inphase across different SS-RTM images. As the number of RTM runs (i.e.,iterations) increases, the probability that an interference inducednoise event does not change phase drops exponentially.

With reference now to FIG. 2, a flow diagram of a method 200 forgenerating an enhanced image is provided in accordance with at least oneembodiment of the present invention. The method 200 may beadvantageously implemented in connection with any appropriate method,such as the method 100, and/or system, such as a computer implementedsystem, to meet the design criteria of a particular application. Themethod 200 generally includes a plurality of blocks or steps that may beperformed serially. As will be appreciated by one of ordinary skill inthe art, the order of the steps shown in FIG. 2 is exemplary and theorder of one or more steps may be modified within the spirit and scopeof the present invention. Additionally, the steps of the method 200 maybe performed in at least one non-serial (or non-sequential) order, andone or more steps may be omitted to meet the design criteria of aparticular application. Block 202 represents an entry point into themethod 200.

At block 204 a plurality of SS-RTM images are generated using different(i.e., unique) random reciprocal sets of encoding functions. In at leastone embodiment the SS-RTM images are generated using the method 100 ofFIG. 1, above. However, any appropriate method and/or system may be usedto generate the images to satisfy the design criteria of a particularapplication.

At block 206 each SS-RTM image is decomposed into its components. In atleast one embodiment, the decomposition is performed by transforming theimage into a different domain (e.g., frequency). In one embodiment alinear transform such as the Fourier transform, wavelet transform,curvelet transform, F-K transform, radon transform, and/or the like maybe used to decompose the image into transform domain coefficients(fourier/wavelet/curvelet, F-K, radon coefficients). However, anyappropriate technique may be used to satisfy the design criteria of aparticular application.

At block 208, image components with substantially similarcharacteristic(s), such as phase and/or magnitude coefficients, acrossthe set of SS-RTM images are identified. Substantially constantcomponents are generally associated with signal energy and, therefore,are retained. Alternatively, if the characteristic(s) of an imagecomponent varies across the set of SS-RTM images, then such a componentpredominantly contains noise energy and, therefore, may preferably beattenuated.

At block 210 an enhanced subsurface image may be generated by combiningthe retained components. In general, the enhanced image will exhibitless noise than an image resulting directly from SS-RTM.

Block 212 represents an exit from the method 200.

In summary, each original SS-RTM image may include constant signal termsas well as variable noise terms. If, for example, a reflection event inan SS-RTM image maintains a substantially constant amplitude acrossSS-RTM images (i.e., runs), then that reflection event is likely to be apart of the desired signal. In contrast, if a reflection event variesacross SS-RTM runs, then that reflection event is a component ofinterference noise and should be attenuated.

In at least one specific embodiment, the method 200 of FIG. 2 mayinclude the following steps:

-   1. A plurality of SS-RTM images are generated using different known    random reciprocal unit-magnitude encoding functions (e.g., step    204).-   2. The curvelet transform of all the SS-RTM images are generated    (e.g., step 206).-   3. The set of curvelet coefficients from all the SS-RTM coefficients    that contain a specific dip, frequency, and location are identified    and combined to form an enhanced curvelet coefficient (e.g., step    208). More specifically:    -   i. The set of curvelet coefficients are averaged to provide a        preliminary signal curvelet coefficient estimate.    -   ii. The noise level in the preliminary signal coefficient        estimate is deduced by computing the variance of the set of        curvelet coefficients    -   iii. If the signal curvelet coefficient estimate's energy is        substantially larger than the noise level, then the coefficient        is retained almost unaltered (primarily signal). If the noise        level is large, then the coefficient is attenuated (primarily        noise). It may be noted that one skilled in the art having the        benefit of this disclosure may select among any one or more of        several techniques to perform the signal to noise comparison        (e.g., hard-thresholding, soft-thresholding, empirical Bayes        thresholding, Wiener filtering, or the like).-   4. The inverse curvelet transform of the retained coefficients may    then be generated to obtain the interference-free RTM image (e.g.,    step 210).

While the present invention has generally been described with encodingperformed using scalars in the time domain, it should be understood thatthe invention may also be extended to include frequency-domainencodings. Different sets of reciprocal source and receiver randomencoding functions, a_(n)(i)[f]'s and b_(n)(i)[f]'s respectively, may beused for each frequency f of the source function and receivermeasurement, such that a_(n)(i)[f]=1/b_(n)*(i)[f] ora_(n)(i)[f]=b_(n)(i)[f]=0. Crosscorrelation in the time domain isequivalent to multiplication with complex conjugation in the frequencydomain. Therefore, by using different reciprocal encodings in thefrequency domain, the SS-RTM image resulting from frequency component fis given by:

$\begin{matrix}\begin{matrix}{{{Image}_{n}\lbrack f\rbrack} = {F\left\lbrack {\sum\limits_{i}{{{a_{n}(i)}\lbrack f\rbrack}{{S(i)}\lbrack f\rbrack}}} \right\rbrack}} \\{{B\left\lbrack {\sum{b_{n}*{(k)\lbrack f\rbrack}R*{\left( {k,j} \right)\lbrack f\rbrack}}} \right\rbrack}} \\{= {\begin{pmatrix}{\sum\limits_{i}{{{a_{n}(i)}\lbrack f\rbrack}b_{n}*}} \\{{(i)\lbrack f\rbrack}{F\left\lbrack {{S(i)}\lbrack f\rbrack} \right\rbrack}{B\left\lbrack {{R\left( {i,j} \right)}*\lbrack f\rbrack} \right\rbrack}}\end{pmatrix} +}} \\{{\begin{pmatrix}{\sum\limits_{i,k}{{{a_{n}(i)}\lbrack f\rbrack}b_{n}*}} \\{{(k)\lbrack f\rbrack}{F\left\lbrack {{S(i)}\lbrack f\rbrack} \right\rbrack}{B\left\lbrack {R*{\left( {k,j} \right)\lbrack f\rbrack}} \right\rbrack}}\end{pmatrix}.}} \\{= {\left( {\sum\limits_{i}{{F\left\lbrack {{S(i)}\lbrack f\rbrack} \right\rbrack}{B\left\lbrack {{R\left( {i,j} \right)}*\lbrack f\rbrack} \right\rbrack}}} \right) +}} \\{{\begin{pmatrix}{\sum\limits_{i,k}{{{a_{n}(i)}\lbrack f\rbrack}b_{n}*}} \\{{(k)\lbrack f\rbrack}{F\left\lbrack {{S(i)}\lbrack f\rbrack} \right\rbrack}{B\left\lbrack {R*{\left( {k,j} \right)\lbrack f\rbrack}} \right\rbrack}}\end{pmatrix}.}} \\{= {{{Signal}\lbrack f\rbrack} + {{Noise}_{n}\lbrack f\rbrack}}}\end{matrix} & \lbrack{EQ12}\rbrack\end{matrix}$

The total image may be generated by summing over all frequencycontributions. Again, as described earlier, unlike the signal component,the noise does generally vary across SS-RTM runs.

A known special case of SS-RTM employs so-called plane waves. In such acase, the source functions and receiver measurements are delayed so thatthe effective source functions and receiver measurements resemble planarwaves. Plane-wave decomposed source functions and receiver measurementsmay then be used as effective sources and receivers and a plane-waveSS-RTM image may be constructed using, for example, equation [EQ12].Note that time-domain delays is equivalent to linear phase encodings inthe frequency domain.

In at least one embodiment, then, step 204 of method 200 may beimplemented using different reciprocal random encoding functions onplane waves with different angle of incidence. Such an approach ismathematically equivalent to constructing the frequency domain encodings(i.e., a_(n)(i)[f]) by randomly combining linear phase encodings. Insuch a case, b_(n)(i)[f]=1/a_(n)(i)[f].

With reference now to FIGS. 3A-3D, images from a specific example areshown to further illustrate the application of one or more embodimentsof the present invention. It should be understood that these images andthe corresponding discussion that follows are provided for clarity andare not intended to limit the present invention in any manner beyond thelimits setout in the attached claims. The individual shots weresimulated with marine streamer geometry.

FIG. 3A shows a simulated stacked RTM image obtained by migrating 383individual shots along the 2D line. As such, this image represents theresults of conventional SS-RTM migration and may be used as a benchmarkto exemplify the quality of the present invention. It may be noted thatthese individual shots were encoded using time-domain binary (±1)encoding scheme.

FIG. 3B shows an image obtained by simulating RTM on one encodedsimultaneous source super-gather. One may note that the image is noisywith poor signal-to-noise ratio (SNR). FIG. 3C shows a correspondingstack of 25 such SS-RTM images, each obtained with a different encoding(generally corresponding to step 204 of method 200). The image is animprovement over the FIG. 3B. However, even after 25 stacks visiblenoise in still observable. Next we apply steps 206 and 208 of method 200using curvelet transform and wiener filtering. The resulting enhancedimage (e.g., result of step 210 of the method 200) is shown in FIG. 3D.The quality of the image is very similar to the image in FIG. 3A; whichis the desired result.

With reference now to FIG. 4A, a distribution of the noise spatially forthe curvelet coefficients of the example of FIGS. 3A-3D is shown. Fromthe figure, one may observe that the noise is greater in the shallowpart than the deeper part. Recognizing such spatial distribution isbeneficial in deciding the threshold used to carry out de-noising. FIG.4B shows the noise distribution in spatial wavenumber domain. Examiningthe distribution in different domains may be beneficial in designingappropriate filters that protect the signal.

As previously mentioned, the present invention has been described in thecontext of SS-RTM. It should be understood, however, that the inventionhas broader applicability and may be used in connection with many othercomputation intensive operations, including iterative inversiontechniques. For example, the present invention may also be applied toefficiently compute gradients during Full Waveform Inversion (i.e.,FWI). In general, FWI estimates a model of subsurface parameters (suchas wave propagation velocities) by iteratively minimizing the differencebetween observed data and data simulated with the subsurface parametermodel. One step in FWI is computing the gradient of a cost function.This gradient computation step is nearly identical to RTM in terms ofbeing computationally expensive. In at least one embodiment the presentinvention may be applied to speed up gradient computation as follows:

-   -   1. Generate the gradients of objective functions computed using        different sets of known random reciprocal encoding functions        (similar to step 204 of method 200).    -   2. Decompose each gradient image into its components by, for        example, transforming each gradient into a different domain        (similar to step 206 of method 200).    -   3. If the characteristics (such as phase and magnitude) of an        image component are similar across all gradients, then the        component contains predominantly signal energy and, therefore,        is retained. In contrast, if the characteristics (such as phase        and magnitude) of an image component vary across all or most        gradients, then the component contains predominantly noise        energy and, therefore, is attenuated (similar to step 208 of        method 200).    -   4. Construct an enhanced gradient by combining all of the        retained components (similar to step 210 of method 200).

Accordingly, one of ordinary skill in the art may appreciate, given thebenefit of this disclosure, that the method 200 may be implemented withany of a number of computationally intensive processes substituted forSS-RTM.

In the simultaneous source method, individual shots are encoded with anencoding function and summed to form a simultaneous shot gather. In caseof fixed-spread geometry, a receiver listens to all shots whereas formarine streamer a receiver listens to only a subset of shots.Unfortunately, in the process of forming these encoded shots, the offsetinformation in the data is lost. This is because a single trace hascontribution from many shots. Therefore the traditional techniques toform image gathers do not work for encoded simultaneous source datasince the source-receiver pair information for the traces is lost. As aresult, offset-based imaging algorithms such as Kirchhoff or Beammigration cannot be applied for simultaneous source data. Despite thisdrawback, it is possible to produce an image using simultaneous sourcedata—since it requires only a forward and backward propagation with acomplicated source function (i.e., the simultaneous encoded source). Thepropagation mechanism can be the 1-way or 2-way wave equation method. Itis also important to note there is no complication if the simultaneoussource data are generated from non-fixed receiver spread such as marinestreamer data since there is no data fitting process in a migrationalgorithm, unlike FWI. Multiple realizations may be used to reducecross-talk noise that occurs in forming these images, with the noisereduction either via simple stacking or using the de-noising techniquedisclosed previously, for example using the curvelet transform.

These principles can be used to develop a velocity model buildingembodiment of the present invention that can be described as follows:

-   (1) Form shot gathers from the data and bin the shot gathers into    offset bins, i.e. shot gathers with a specified range of offsets.    Possible offset bins can be near, mid and far but any possible    combinations are permissible.-   (2) Encode each offset bin to form simultaneous source data.-   (3) Use M realizations of a particular offset-binned simultaneous    source data to generate the image. This can be achieved by RTM or a    1-Way wave propagation method. Different realizations means using a    different random set of encoding functions. Each realization is    migrated and stacked, thereby improving the signal-to-noise ratio.    Alternatively, the de-noising technique disclosed previously, for    example using the curvelet transform, could be used on the multiple    realizations.-   (4) If the velocity model is accurate, one would expect that events    in the individual images from different offsets will be registered    at the same depth. Any mis-positioning of these images indicates    inaccuracy in the velocity model. This mis-positioning information    may be translated into a velocity update using any of several    approaches, for example:    -   (a) Determine the regions where the images are mis-positioned;        these will be the regions where local velocity updates are        needed to properly position the events in the offset binned        migrated images. Different velocity panels with local        modifications can be used to regenerate the images and examine        the image registration.    -   (b) Another possibility is that images from mid and far offsets        can be cross-correlated to a near-offset image to determine the        shift. These shifts then can to be translated into velocity        update.    -   (c) Treat the image mis-positioning as 3D image registration        problem. In this case, offset gathers can be produced and the        moveout can be used to update the velocity using reflection        tomography.        The cost of this approach is N×M where N is the number of offset        bins (typically ˜3) and M is the number of realizations used for        each offset bin simultaneous source imaging. In regions with        large number of shots, this can be computationally efficient.        The SS-RTM images can also be used to compute angle gathers        using either non-zero lag cross correlations or using Poynting        vectors. A key idea is that instead of migrating individual        shots, one can use simultaneous source RTM to compute the angle        gathers.

FIG. 5 shows the RTM images of synthetic data binned into differentoffset bins. The RTM was carried out with sequential shots, which is 20times more expensive than the SS-RTM approach shown in FIG. 6. The firstpanel in FIGS. 5 and 6 is the RTM image for offsets from 0-1000 m; thesecond panel is the RTM image for offsets from 1500-2500 m; the thirdpanel is the RTM image for offsets from 2000-3000 m; and the fourthpanel is the RTM image for offsets from 2500-3500 m.

FIG. 6 shows similar RTM images, but using SS-RTM. The close equivalenceof FIGS. 5 and 6 illustrates that SS-RTM is able to produce verycomparable images compared to sequential RTM at much reduced cost. Thespeed up generated by simultaneously migrating several sourcessimultaneously using source encoding was taught by Ober and Romero (U.S.Pat. No. 6,021,094 to Ober, et al.). However, the present inventivemethod adds a denoising strategy to further speed up the migration andalso adds the idea of interrogating the velocity model using offsetbinned migrated images. One can notice from both FIG. 5 and FIG. 6 thatwith increasing depth, focusing of the images is lost as offsetsincrease, and careful examination shows that the some of the deeperevents are mis-positioned indicating that velocity in some of the deeperregions can be incorrect.

FIG. 7 shows an extreme example where error is deliberately introducedin the velocity model. The RTM images for the various offset bins withthe incorrect velocity model clearly indicate the mis-positioning aswell as loss in signal strength (i.e. loss of focusing) with increasingoffset.

In accordance with various embodiments of the present invention, themethods described herein are intended for operation as software programsrunning on a computer processor. Dedicated hardware implementationsincluding, but not limited to, application specific integrated circuits,programmable logic arrays, microprocessor based devices, and otherhardware devices can likewise be constructed to implement the methodsdescribed herein. Furthermore, alternative software implementationsincluding, but not limited to, distributed processing orcomponent/object distributed processing, parallel processing, or virtualmachine processing can also be constructed to implement the methodsdescribed herein.

It should also be noted that the software implementations of the presentinvention as described herein are optionally stored on a tangiblestorage medium, such as: a magnetic medium such as a disk or tape; amagneto-optical or optical medium such as a disk; or a solid statemedium such as a memory card or other package that houses one or moreread-only (non-volatile) memories, random access memories, or otherre-writable (volatile) memories. A digital file attachment to email orother self-contained information archive or set of archives isconsidered a distribution medium equivalent to a tangible storagemedium. Accordingly, the invention is considered to include a tangiblestorage medium or distribution medium, as listed herein and includingart-recognized equivalents and successor media, in which the softwareimplementations herein are stored.

While embodiments of the invention have been illustrated and described,it is not intended that these embodiments illustrate and describe allpossible forms of the invention. Rather, the words used in thespecification are words of description rather than limitation, and it isunderstood that various changes may be made without departing from thespirit and scope of the invention.

REFERENCES

-   Baysal, E., Kosloff, D. D., and Sherwood, J. W. C., 1983, Reverse    time migration, Geophysics, 48, 11, 1514-1524.-   Dai, W., and Schuster, J., 2009, Least-squares migration of    simultaneous sources data with a deblurring filter, SEG expanded    abstract, 2990-2994.-   Godwin, J, and Sava, P., 2010, Simultaneous source imaging by    amplitude encoding, CWP report 645.-   C. C. Ober, L. A. Romero, and D. C. Ghiglia, Method of migrating    seismic records, U.S. Pat. No. 6,021,094 (assignee: Sandia)-   L. A. Romero, D. C. Ghiglia, C. C. Ober, and S. A. Morton, 2000,    Phase encoding of shot records in prestack migration, Geophysics,    vol. 65, 426-436.-   F. Perrone and P. Sava, Comparison of shot encoding functions for    reverse-time migration, 2009, SEG Expanded Abstract 2980-2984.-   F. Perrone and P. Sava, 2010, Wave-equation migration with dithered    plane waves, CWP report 646.-   Claerbout, J. F., 1985, Imaging the earth's interior: Blackwell    Scientific Publications.-   Stolt, R. H., 1978, Migration by Fourier transform: Geophysics, 43,    23-48.-   Tang, Y., and Biondi, B., 2009, Least-squares migration/inversion o    blended data, SEG Expanded abstract, 1041-1044.-   Verschuur, D. J., and Berkhout, A. J., 2009, Target oriented    least-squares imaging of blended data, SEG expanded abstract,    2889-2893.-   Whitmore, N. D., 1983, Iterative depth migration by backward time    propagation, SEG Expanded Abstract, 382-385-   Yilmaz, 1987, Seismic Data Processing, Society of Exploration    Geophysicists.

1. A method for processing seismic data, the method comprising: a.obtaining a plurality of initial subsurface images, wherein each of theinitial subsurface images is generated using a unique random set ofencoding functions; b. decomposing each of the initial subsurface imagesinto components; c. identifying a set of components comprising one of(i) components having at least one substantially similar characteristicacross the plurality of initial subsurface images, and (ii) componentshaving substantially dissimilar characteristics across the plurality ofinitial subsurface images; and d. generating an enhanced subsurfaceimage using the set of components identified in step c; wherein at leastone of steps a-d is performed using a computer.
 2. The method of claim 1wherein the initial subsurface images and the enhanced subsurface imageare SS-RTM images.
 3. The method of claim 1 further comprising the stepof attenuating components not identified as having at least onesubstantially similar characteristic across the plurality of initialsubsurface images.
 4. The method of claim 1 wherein decomposing eachinitial subsurface image includes transforming each initial subsurfaceimage into another domain.
 5. The method of claim 4 wherein the domainis frequency.
 6. The method of claim 4 wherein a curvelet transform isused to transform each initial subsurface image into another domain. 7.The method of claim 4 wherein at least one of a Fourier transform,wavelet transform, F-K transform, curvelet transform, and radontransform is used to transform a subsurface image into another domain.8. The method of claim 1 wherein the initial subsurface images aregenerated by: a. obtaining a set of shot gathers comprising forward andbackward wave component data; b. selecting first and second randomencoding functions; c. encoding the forward wave component data for eachsource in the set of shot gathers using the first random encodingfunction to form an Encoded Source Super-Shot Wave Component; d.encoding the backward wave component data for each receiver in the setof shot gathers using the second random encoding function to form anEncoded Receiver Super-Shot Wave Component; e. forward propagating theEncoded Source Super-Shot Wave Component to generate a ForwardPropagated Wave Component; f. back propagating the Encoded ReceiverSuper-Shot Wave Component to generate a Back Propagated Wave Component;g. applying an imaging condition to the Forward and Back Propagated WaveComponents to generate a subsurface image; and h. iteratively repeatingsteps b-g until a predetermined condition is satisfied, wherein thefirst and second random encoding functions are selected so that thefunctions are unique for each iteration.
 9. The method of claim 8wherein the encoding the forward and backward wave components isperformed using scalars in the time domain.
 10. The method of claim 8wherein the encoding the forward and backward wave components isperformed using scalars in the frequency domain.
 11. The method of claim8 wherein the first and second random encoding functions are reciprocal.12. The method of claim 11 wherein at least one of the first and secondreciprocal random encoding functions is a unit-magnitude encodingfunction.
 13. The method of claim 12 wherein one or more of the firstand second reciprocal random encoding function include a unit-magnitudecomplex number encoding function.
 14. The method of claim 8 wherein thefirst and second random encoding functions include reciprocal randomencoding functions on plane waves with different angles of incidence.15. The method of claim 8 wherein the first random encoding function isequivalent to the second random encoding function.
 16. A method forprocessing seismic data, the method comprising: a. obtaining a pluralityof SS-RTM subsurface images, wherein each of the SS-RTM subsurfaceimages is generated using a unique random unit-magnitude set of encodingfunctions; b. decomposing each of the plurality of SS-RTM subsurfaceimages into curvelet coefficients; c. averaging the curveletcoefficients to generate a preliminary signal curvelet coefficientestimate; d. computing a variance of a subset of the curveletcoefficients to determine a noise level in the preliminary signalcurvelet coefficient estimate; e. attenuating noise in the curveletcoefficients using the determined noise level and the preliminary signalcurvelet coefficient estimate to generate attenuated curveletcoefficients; and f. performing an inverse curvelet transform on theattenuated curvelet coefficients to generate an enhanced SS-RTM image;wherein at least one of steps a-f is performed using a computer.
 17. Themethod of claim 16 wherein step e is performed using one or more ofhard-thresholding, soft-thresholding, empirical Bayes thresholding, andWiener filtering.
 18. A method for computing a gradient of a costfunction associated with seismic data, the method comprising: a.generating a plurality of gradients of objective functions computedusing unique sets of random reciprocal encoding functions; b.decomposing each of the gradients into components; c. identifying a setof components having at least one substantially similar characteristicacross the plurality of gradients; and d. generating an enhancedgradient using the set of components identified in step c; wherein atleast one of steps a-d is performed using a computer.
 19. The method ofclaim 1, wherein the initial subsurface images are obtained bysimultaneous-source reverse-time migration, the encoding functions haveunit magnitude, the decomposing is performed by curvelet transformgenerating a set of curvelet coefficients for each initial subsurfaceimage, and step c comprises: averaging the curvelet coefficients togenerate a preliminary signal curvelet coefficient estimate; andcomputing a variance of at least a subset of the curvelet coefficientsto determine a noise level in the preliminary signal curveletcoefficient estimate; and step d comprises: attenuating noise in thecurvelet coefficients using the determined noise level and thepreliminary signal curvelet coefficient estimate to generate attenuatedcurvelet coefficients; and further comprising performing an inversecurvelet transform on the attenuated curvelet coefficients to generatean enhanced simultaneous-source reverse-time migration image.
 20. Themethod of claim 1, further comprising first iteratively inverting theseismic data, said inversion involving computing gradients of objectivefunctions associated with the seismic data, and then performing themethod with the gradients being regarded as the initial subsurfaceimages.
 21. The method of claim 1, further comprising initially formingshot gathers from the seismic data and binning the shot gathers into atleast two bins, each offset bin having shot gathers with a specifiedrange of offsets; wherein step a comprises: encoding the gathers in eachoffset bin using a unique random set of encoding functions to formcomposite gathers of simultaneous source data, then repeating at leastonce using a different random set of encoding functions, thereby formingat least two realizations of each offset-bin composite gather, said atleast two realizations becoming, after migration using an assumedvelocity model, the plurality of initial subsurface images for step b.22. The method of claim 21, wherein the migration is SS-RTM, and furthercomprising examining coherency or consistency of the enhanced subsurfaceimage for different offset bins to assess accuracy of the assumedvelocity model.
 23. The method of claim 22, further comprising usingmis-positioning of one or more reflection events between differentoffset bins to estimate a corresponding update to the assumed velocitymodel.
 24. A method for processing seismic data comprising: forming shotgathers from the data and binning the shot gathers into at least twobins, each offset bin having shot gathers with a specified range ofoffsets; encoding the gathers in each offset bin using a unique randomset of encoding functions to form composite gathers of simultaneoussource data, then repeating at least once using a different random setof encoding functions, thereby forming at least two realizations of eachoffset-bin composite gather; migrating and stacking the at least tworealizations of each offset-bin composite gather, wherein an assumedvelocity model is used to generate seismic images using reverse timemigration or a one-way wave propagation method; and assessing accuracyof the assumed velocity model based on whether a reflection event ispositioned at a same depth in different offset-bin images; wherein atleast the encoding and the migrating and stacking are performed using acomputer.
 25. The method of claim 24, wherein the migration is SS-RTM,and further comprising translating mis-positioning of reflection eventsinto an update to the assumed velocity model.